Thin superconductors and SQUIDs in perpendicular magnetic field 
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It is shown how the static and dynamic electromagnetic properties can be calculated for thin flat 
superconducting films of any shape and size, also multiply connected as used for SQUIDs, and for 
any value of the effective magnetic London penetration depth A. As examples, the distributions of 
sheet current and magnetic field are obtained for rectangular and circular films without and with 
slits and holes, in response to an applied perpendicular magnetic field and to magnetic vortices 
moving in the film. The self energy and interaction of vortices with each other and with an applied 
magnetic field and/or transport current are given. Due to the long ranging magnetic stray field, 
these energies depend on the size and shape of the film and on the vortex position even in large 
films, in contrast to the situation in large bulk superconductors. The focussing of magnetic fiux into 
the central hole of square films without and with a radial slit is compared. 

PACS numbers: 74.78.-w, 74.25.Ha,74.25.0p 



I. INTRODUCTION 

The calculation of the electromagnetic properties 
of thin superconducting films of finite size as used, 
e.g., in Superconducting Quantum Interference Devices 
(SQUIDs)i is a complicated problem since these sensi- 
tively depend on the shape and size of the film. This 
is so since the currents in thin films are not screened as 
in bulk superconductors, but interact via the magnetic 
stray field they generate outside the film. In particular, 
the self energy of a magnetic vortex in thin films (called 
Pearl vortey^) and the interaction between two such vor- 
tices depend on their position in the film even for very 
large films. This means a vortex is never "far from the 
film edges" , in contrast to the behavior of vortices in the 
bulk, whose energy, current density, and magnetic field 
become independent of the vortex position and of the 
specimen shape when the distance from the surface is 
much larger than the London penetration depth A. This 
is so since in the bulk the factor 1/(1 -I- A:^A^) in the 
Fourier transforms (with wave vector components k^, ky, 
kz, k\ — k"^ -\- ky) causes the fields and currents to de- 
cay exponentially over the length A at large distances 
from the vortex core. In thin films of thickness d < X 
the effective magnetic penetration depth A = X^/d is 
larger than A and the Fourier transforms contain a factor 
^/{^k± -\- k'^A) that describes also the long-range non- 
exponential interaction of vortices via the magnetic stray 
field outside the film. 

But even films in the ideal Meissner state containing no 
vortices present a difficult problem. Properties of macro- 
scopic circular disks and rings in a perpendicular applied 
magnetic field were calculated recently for ideal screening 
(A = 0)^ and for arbitrary A4 When this ring has a radial 
slit, e.g. in a washer-shaped SQUID, the circular sym- 
metry is lost, but some properties like the sheet current 
and the concentration of magnetic flux into the central 
hole (flux focussing) can still be calculated approximately 
from this circular symmetric model by forcing the current 
in the ring to be zero.^i'* (This situation may be achieved 



by appropriate magnetic history.) Below we shall com- 
pare this approximation with the exact two-dimensional 
(2D) solution for a slitted ring and flnd partial agreement 
(Sec. III). 

While the slitted ring or slitted square fllm with an 
applied magnetic held Ha and/or transport current la 
are simply connected geometries (Fig. 1, right two plots), 
a closed ring or a slitted film with the entrance of the 
slit short-cut by superconducting contacts (e.g. by weak 
links) present multiply connected geometries (Fig. 1, left 
two plots). These are more difficult to calculate since the 
(quantized) magnetic flux <& (or fluxoid <I>f when A > 0) 
trapped in this hole, and the current / circulating around 
the hole, are additional parameters, besides Ha and Ia- 
in films with n holes or slots that are fully surrounded by 
superconducting material, there are n such fluxoids $fi 
and currents li that depend on the magnetic history and 
that may be used to deflne n self-inductances Li — ^a/Ii 
and n{n — l)/2 mutual inductances Mij — ^u/Ij. 

This paper shows how all these (actually 3D) thin-film 
problems can be solved numerically by a 2D matrix in- 
version method allowing for non-equidistant grid points. 
The presented general equations and concrete examples 
generalize previous methods that either work only for 
equidistant spatial grids^ (which are not very accurate 
near the film edges or near a narrow slit or small hole), 
or for general grid did not account for finite penetration 
depth A > 0)^ or assumed simply connected geometry;^ 
or applied only to circular disks or ringsj2i4 In this paper 
I consider the electrodynamics of finite-sized macroscopic 
films that can be described by London theory, which ap- 
plies when the magnetic field is much smaller than the 
upper critical field Hc2 and the superconductor is much 
larger than the coherence length ^. The application to 
SQUIDs will be dealt with elsewhere.-' In the present pa- 
per I shall thus not need the notions of fiuxoid quan- 
tization, phase of the superconducting order parameter, 
vector potential, and voltages caused by the Josephson 
effect, but some quantities computed here will be needed 
in the theory of SQUIDs, e.g., self-inductance and effec- 
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tive area. 

Our 2D matrix inversion method can be quite accurate 
even when the number of grid points is not very large. 
For example, in the ideal Meissner state with A = 0, the 
computed current density exactly yields Hz{x,y,0) — 
inside the film (z = is the film plane). For A > 0, accu- 
rate results are obtained even when A is smaller than the 
spacing of our rectangular grid. For similar calculations 
using a finite-element method and a different integral ker- 
nel (or matrix) see Ref. Isl Analytical^ and numericaliS 
London calculations were performed in the limit of large 
A, i.e., for small disks with radius R A. (Here I shall 
not list numerous recent work on mcsoscopically small 
superconductors with vortices computed from Ginzburg- 
Landau theory.) An elegant and fast method that com- 
putes the currents in films from the magnetic field pat- 
tern measured at the film surface with high resolution, 
without having to store or explicitly invert a large ma- 
trix, is described by Wijngaarden et ai.-fLiJ^ this method 
has all advantages of the direct matrix inversion method 
and avoids the inversion by Fourier transform, that would 
require knowledge of the magnetic field pattern also out- 
side the film area. The static Bean model for thin films of 
any shape is computed by Prigozhin using a variational 
methodii^ 



4. The integral of g{x, y) over the film area equals the 
magnetic moment of the film if g = on its edge. 

5. The difference g{vi) — g{v2) is the current that 
crosses any line connecting the two points ri and r2. 

6. If the film contains an isolated hole or slot such that 
magnetic flux can be trapped in it or a current / can 
circulate around it, then in this hole one has g{x, y) = 
const = / if g{x, y) = is chosen outside the film. 

7. In a multiply connected film with n holes, n inde- 
pendent constants gi . . . g-n can be chosen for the values 
of g{x,y) in each of these holes. The current flowing 
between hole 1 and hole 2 is then gi — 52 • 

8. A vortex with flux $0 in the film moves in the poten- 
tial V = —^og{x, y), since the Lorentz force on a vortex 
is -J X z$o = -«f>oZ X (z X Vg) = $oVg(a;, y) = -W. 

9. A vortex moving from the edge of the film into a 
hole connected to the outside by a slit, at each position 
{x,y) couples a fluxoid g{x,y)^o/I into this hole, where 
g{x,y) is the solution that has g{x,y) = / in this hole 
(with closed slit, see point 6) and g — outside the film. 

10. When the film has n holes, a vortex in the film at 
{x,y) couples a fluxoid gi{x,y)^Q/I into the ith hole (if 
this is connected to the edge by a slit), where gi{x, y) is 
the solution exhibiting gi — I in this hole and gi = in 
all other (isolated) holes and on the fllm edge. 



II. CALCULATION METHOD 

This section describes how for thin flat superconduct- 
ing films of any shape, also multiply connected as needed 
for SQUIDs, one can calculate the static and dynamic 
response to an applied magnetic field, applied electric 
current, and to vortices moving in the film. In such 
problems the central physical quantity is the thickness- 
integrated current density, called sheet current J(x,?/) = 
Jdzj{x,y,z) — {Jx, Jy). For films with constant thick- 
ness d and nearly z-independent current density j(x, y, z) 
one approximately has J = jd, but the following equa- 
tions are more general, applying also to films with spa- 
tially varying thickness (i(x, y) if the typical length of this 
variation is 3> d. 



A. Properties of the stream function 

Since one has V • J = in the film except at small 
contacts, one can express J in terms of a scalar potential 
or stream function g{x, y) as 3 — — z x Vg — V x (zg) — 
{dg/dy, —dg/dx). The function g{x, y) has several useful 
properties and interpretations: 

1. g{x,y) is the local magnetization or density of tiny 
current loops. 

2. The contour lines of g(x, y) are the current stream 
lines. Typical g{x,y) look like a mountain (Fig. 2). 

3. On the boundary of the film one may put g{x, y) = 
since the boundary coincides with a stream line. 



B. Ampere's law for thin films 

From Ampere's local 3D law j = V x H one obtains 
for a current-carrying film in the plane z = a nonlo- 
cal 2D relation between the perpendicular magnetic field 
Hz{x, y, 0) and the stream function g{x, y): 

HAr) = Ha{r)+ f dV g(r, r') 5(r') • (1) 
Js 

Here r = (x,y), Ha{v) is the z component of the applied 
magnetic field, and S is the area of the film. The inte- 
gral kernel Q(r, r') has the meaning of the magnetic field 
(along z) caused at point r in the plane of the film by 
a magnetic dipole (or tiny current loop) of unit strength 
positioned at r' and directed along z. From the known 
dipole field in any plane z = const, one obtains formally: 

Q(r, r') = Q{p) = lim , f f ^ C/2 ' (2) 
z^o 47r(z^ + p'^)°/^ 

Note that this kernel depends only on the distance p — 
|r - r'|. For p 7^ one has Q{p) = -l/iinp^), but the 
integral of Q{p) over the infinite plane vanishes, i.e., the 
total magnetic flux of a dipole is zero in any plane z = 
const, also for z — s- 0, thus Q{p) is highly singular at 
p = 0. For explicit calculations one has to decide how to 
deal with this singularity of Q. 

For numerics, one has to write the integral as a sum. 
Defining a 2D grid whose N points fill the film area 5*, 
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one may approximate any integral over S* by a sum, 



cPrfir) 



N 

E 

1=1 



(3) 



where the Wi are weights, e.g., Wi = const — S/N 
for equidistant grids that avoid the film edges, see Ap- 
pendix A. The accuracy of this numerical integration 
can be strongly increased by choosing an appropriate 
non-equidistant grid, e.g., a grid that is denser near the 
boundary of S or near possible poles or jumps of the 
integrand f{x,y). Equation (1) now becomes 



H^{ri) = Ha{ri) + ^ Wj g{rj) 



(4) 



with the matrix Qij = Q(ri, fj). Equation (4) is formally 
solved by matrix inversion, i.e., by writing 



gin) 
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where Kij is the inverse matrix 



{Q^. 



(5) 



(6) 



defined by the equation Kn^QijWj) = Sij with Sij — 1 
for i = j and Sij — for i ^ j. Note that, in contrast 
to Qij, the inverse matrix Kij depends on the shape of 
the film and not only on the difference — . For the 
film shapes we have tested, all the matrix elements Kij 
are found to be negative, with a sharp negative peak at 
i = j, and tending to zero when or Vj approaches the 
film edge, see Fig. 3 for an example. In principle, the in- 
verse kernel K{x,y;x' ,y'), integrated over y and y' (the 
strip width) was introducedi^ and depictedi^ earlier in 
the context of the magnetostatic energy of a tilted and 
curved narrow superconducting strip with pinned vor- 
tices. Then, and later K{x, y; x' , y') was calculated by 
iterating an integral equation. 

A useful expression for the matrix Qij was obtained in 
terms of a Fourier series in Ref.H but this method works 
only for equidistant grids, while for non-equidistant grids 
the matrix inversion is singular. This numerical form of 
Qij is thus not very accurate when one is interested in 
the sharply peaked J and Hz near the film edges, or 
when the film exhibits fine structures, e.g., a small hole 
or narrow slit, since the maximum number of grid points 
on present Personal Computers is limited to about N — 
5000, yielding a very large N x N matrix Qij. 

A matrix Qij that works well for any grid , also with 
non-constant weights Wi, is obtained as follows. From 
Eq. (2) one has for i ^ j: 



1 



(7) 



The diagonal terms Qn are obtained from the condition 
that the integral of Q{p) taken over the infinite area has 



to vanish. Splitting this integral into the integral over 
the film area S plus the integral over the infinite area S 
outside 5, and writing the first integral as a sum, we get: 

= / (fr'Q{v,-v') 

J OQ 

]Q,jWj+J_d\'Q{r,-r'). (8) 



i 



Defining a 2D function C(r) as an integral over the film 
area or as a contour integral over the film edge, 



Cir) 



dPr' 



Att\y — r' 



d4) 



47ri?(0) 



(9) 



with i? = |r — r'l ((/) is the angle between the vector 
R = r — r' pointing to the point r' on the film edge 
and any fixed direction, say, the -I- a; axis), and writing 
Ci = C{ri), q^j = l/(47r|ri - r^p), we get from (8) the 
diagonal term QaWi — Tlij^i ^ij^j + ^i- The full matrix 
reads thus: 

Qij = (% - 1) Qij + 5ij ^ qu wi + Ci^/wj . (10) 

Note that the terms in (10) should not be rearranged 
since qa — oo. The matrix (10) is well behaved during 
inversion, so one may write 



K- 



((5y - 1) qijWj + 5ij ^ qu wi + Q 



(11) 



For a rectangular film filling the area \x\ < a, \y\ < b 
one has explicitly from Eq. (9)»^^ 



Cix,y) 



(b ~ qvY 



1/2 



(12) 



with p,q — ±1. Interestingly, expression (12) may be 
used also for films that have a hole or slit, or several 
holes, or that do not fill the rectangle completely, e.g., 
a circular disk with radius < a ^ b. In such cases one 
has to omit in Eq. (4) and (5) the grid points that fall 
outside the film (but keep the points in isolated holes, 
see Set. II E). Therefore, Eq. (5) with the explicit kernel 
Kij from Eq. (11) and Ci from Eq. (12), allows one to 
compute the stream function g{x,y) and thus the sheet 
current for thin films of arbitrary shape when one knows 
the magnetic field Hz(x,y) — Ha{x,y) generated inside 
the film by this current. Information on outside the 
film is not required for this stable inversion method. 



C. Static solution of London equation 

In superconductor films with thickness d < A, the Lon- 
don penetration depth, the 3D static London equation 
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A^Vxj + H = may be integrated over the film thick- 
ness d to yield the 2D equation 

H^Xx, y) = -A[V X J{x, y) ] z = AV^gix, y) , (13) 

where A = A^/d is the effective London depth of the film. 
Eliminating Hz{x,y) from Eqs. (1) and (13) one obtains 
an implicit equation for g{x,y): 

Ha{r) = -( dPr' Q{r, r') g(r') + AV^g{r) 
Js 

= -( d?r' [g(r, r') - 5[v - y')KV^] g(r') (14) 
Js 

or with the discretized Eq. (4), 

Hair.) = - E ( - ^"^l ) ■ (15) 

j 

In it the matrix Vf^ computes the Laplacean = 
jdx^ + /dy^ at r = of a function defined on a 
grid, e.g., from the values g{vj) at = and its four 
nearest neighbors. Equation (15) is solved for g(x,y) by 
matrix inversion: 

g{v,)^-Y,Kt,Ha{y,) (16) 

3 

with the inverse matrix 

now depending on A. This matrix inversion is the more 
stable the larger is A, since finite A increases the diagonal 
terms and makes the resulting K^^ a smoother function 
as compared to the case A = considered in Eq. (6). Ex- 
amples for Hz{x, y) are shown in Fig. 4 for a square with 
slit and hole, while Fig. 5 shows some profiles Hz{XtG) 
along the y-axis for the same square. Note that even 
small A = O.Ola allows to partly penetrate the entire 
film. 



D. Dynamic solution of London equation 

The time dependent behavior of superconducting films 
containing vortices may be described within continuum 
theory by the following realistic relation between the lo- 
cal electric field E(a;,y,t) and the sheet current J and 
magnetic induction B — /ioH:'^ 

E = p,(J,B)J(r,t)+/ioAj(r,t). (18) 

Here ps = p/d is the sheet resistivity caused by mov- 
ing vortices, and the second term with A = A^/d and 
J = 9 J/9i is the London term describing acceleration of 
Cooper pairs. The isotropic model (18) assumes that the 
resistivity p depends only on the magnitudes J and B. 
For example, without vortex pinning and Hall effect, one 
has free flux flow with p = ppF ~ {B / Bc2)Pn, where p„ 



is the resistivity in the normal conducting state and Bc2 
is the upper critical field. For thermally activated de- 
pinning a realistic model is p — pq\J I Jc{B)\'^ with creep 
exponent a ^ \ and an in general B dependent critical 
sheet current Jc{B). For a generalization to anisotropic 
superconductors see Ref. |3 

From the induction law B = —V x E, which in the film 
plane reduces to B^ — dE^/dy — dEy/dx, and from J = 
— z X Vg, one obtains poHz = Bz ~ V[psV.g] + poA\7'^g. 
Inserting this into the time derivative of Eq. (1) one finds 
an equation for g{r, t): 

f d^r'[Q{r,r') - 6{r - r')AV^]g{r' ,t) 
Js 

= f{r,t)^Ha{r,t), 
f{v,t) = p,'W[ps{r,t)Wg{r,t)] . (19) 
In discretized form this becomes [cf. Eq. (15)]: 

E(Qy - AV?^.) g{r,) = f{r,,t) - Ha{r^,t) . (20) 

j 

Inverting this one obtains the equation of motion for 
g{x,y,t) in explicit form: 

5(r„ t) = ^ ^ [ fiy, ,t)-Ha (r, , t) ] (21) 

j 

with K:lj from Eq. (17). In the Meissner state or for 
rigidly pinned vortices one has ps — and these dynamic 
equations reduce to the static equations of Set. II C. 

E. Multiply connected films 

Multiply connected films have one or more holes or 
slots that are completely surrounded by superconducting 
material and thus can trap magnetic flux. As a fun- 
damental example I consider here a fllm containing one 
such hole with flux trapped such that a current / circu- 
lates around the hole when no magnetic field is applied. 
Ha — 0. In this case one has g{x, y) = outside the film 
and g{x,y) = / in the hole, and inside the film g{x,y) 
smoothly goes from to I. Generalization of this ex- 
ample to the presence of more holes and to Ha 7^ is 
possible by linear superposition. 

This problem may be solved in three steps. First, con- 
sider the situation where g = I m. the hole and (7 = 
everywhere outside the hole. This means a sharply lo- 
calized sheet current of size / flows along the edge of the 
hole where this g{x, y) has a jump. Such an edge current 
formally can be caused by an effective applied field 

Hf{v,)^-I J2 (ft.^'.-AVy, (22) 

j in hole 

cf. Eq. (15). Next the real sheet current in the film is 
found as the J = — z x Vg that generates this iJ^*(r) 
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inside the film, cf. Eq. (16): 



slit closed, H = 0, I > 0, trapped flux slit open, H > 0, 1 = 0, flux focussing 



5(r) 



- ^ Hf{rj) for Vi in the film, 

j in film 

I for Vi in the hole, 







outside the film. (23) 



Finally, the real magnetic field in the entire plane z = 
is obtained as the field caused by this current, cf. Eq. (4): 



E 

3 



QijWjg{rj), 



(24) 



where now the sum is over all rj in the film and in the 
hole. This method works for any value of A and yields 
continuous functions g{x,y) and Hz{x,y) when A > 0. 
When A = (ideal screening) the resulting Hz{x,y) 
has sharp jumps at the film edges since in that case 
Hz exactly vanishes inside the film and has sharp infini- 
ties just outside the film edges (also in the hole), where 
Hz oc {6 is the distance from the edge). The current 
density in the film for A = diverges similarly, but on the 
inner side of the edges, where g oc (5^/^ and J oc S~^^'^. 
For A > the sheet current at the edges is finite and the 
infinities of H^ are logarithmic on both sides of all edges. 



F. Individual vortices in the film 

The sheet current, magnetic field, and energy of vor- 
tices in the film is obtained by linear superposition from 
the solution for one vortex and the interaction energy of 
a vortex pair. In a film of finite size all these results ex- 
plicitly depend on the vortex positions and not only on 
their distances due to the strong effect of the film edges. 
The existence of one vortex at position r„ modifies the 
static London Equation of Set. II C to give A^V x j + H = 
($o/A<o)z'^2(r-r„), where $o = h/2e = 2.07- IQ-^^ Tm^ 
is the quantum of flux and S2 (r) is the 2D delta function. 
Equations (13) and (14) then become 

JT,(r) - AV^g{r) = Mo S2{r - , (25) 

/ cfr'[Q{v,r')-6{v-r')AV^]g{r') 
Js 

= 110^^0 52{r-r,)-Ha{r). (26) 
Writing the integral as a sum yields 

^(g,,«;, -AVy5(r,) 
j 

= 110^^0 52{ri-r,)-Ha{ri). (27) 

To invert this and find g{ri) and the vortex interaction, 
we have to assume that the vortex sits on a grid point, 
r„ = Tj. Averaging over the grid cell centered at rj 
and having an area Wj, replaces 62(^,1 — Vj) by Sij/wj. 
Inverting this and performing the sum containing Sij then 




FIG. 1: Current stream lines in a thin film square with square 
hole and radial slit, and in a circulax disk with circular hole 
and slit, in the ideal Meissner state A = 0. Top left: Slit 
bridged at the edge, circulating current / > flows due to 
flux trapped in the hole and slit, no applied field Ha = 0. 
Top right and bottom right: Slit open, applied field Ha > 0, 
magnetic flux enters the slit and is focussed into the hole 
where H{x,y) > Ha- Bottom left: Closed slit, applied field 
Ha > 0, some current / > flows such that the flux in hole 
and slit is exactly zero (ideal screening); this state is a su- 
perposition of the two upper states. Note that the current 
near the hole circulates in opposite direction, except in the 
trapped-flux case. 



yields the stream function caused by a vortex positioned 
at Tj and by the apphed field Ha{ri): 



gin) = Mo^^oi^^^M -^K^H^in) 



(28) 



with the inverse matrix Kij from Eq. (17), see also 
Eq. (16). 

A second vortex sitting at sees the potential 
— $oi?(ri); the interaction energy between two vortices 
positioned at rj and is thus: 



(29) 



This potential is repulsive (positive) and sharply peaked 
at Tj = Yj, since all the K^- are negative. One can 

show that the matrix K^/wj is indeed symmetric in 



K^Jw 



iji -^3 — ^fih'^i- For A = this is directly seen from 



the definition, Eq. (6), writing 5ij = J2i ^niQij^j) — 
'^i{Kii/wi){wiQijWj), which shows that Ku/wi is the 
inverse of the symmetric matrix WiQijWj and is thus sym- 
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1 -1 



FIG. 2: The three examples of Fig. 1 for the stream function 

g{x, y) in a square fikn (size 2 x 2 in units of half width a) with 
square hole (half width ai = 0.2a) and open slit (width 0.04 
a), for 40 x 40 grid points, as Fig. 3. For penetration depth 
A = (ideal screening). At the film edges g{x,y) goes to zero 
with vertical slope, and outside the film g = Q (for A > 0, the 
slope of g at the edges is finite, equal to the sheet current). 
Top: Constant applied field Ha ~ 1, open slit meaning a 
current / = 0, thus g — in slit and hole (like Fig. 1 top 
right). Middle: Ha = 0, current 1=1 flowing around the 
hole and closed slit due to trapped flux, yielding p = 1 in hole 
and slit (like Fig. 1 top left). Bottom: Top and middle cases 
superimposed such that the flux trapped in the hole and slit 
is zero (like Fig. 1 bottom left): weights 1 and 1.2, thus in slit 
and hole g = I = 1.2aHa. 




FIG. 3: Example for the inverse matrix K^j = K^{ri,rj) 
in a square film (half width a) with square hole (half width 
oi = 0.2a) and open slit (width 0.04a), similar to the squares 
in Figs. 1 and Fig. 2. For 40 x 40 grid points r; and constant 
Tj = (0.48, 0.42)a. penetration depth A = 0.1a. At the film 
edges and in hole and slit one has K{) = 0. The plotted 
is approximately syunnetric in i. j since here the weights 
Wi w const. It also shows the interaction between a vortex 
at {x,y) and a vortex pair (due to the imposed symmetry) 
sitting at {xj,yj) = (0.48, ±0.42)o, cf. Eq. (29). Its contours 
are the current stream lines of this vortex pair. For A = 0, 
Kf"j = Kij looks similar, but the peak is higher and the slopes 
at the edges are vertical. 

metric itself. For A > one also has to prove that the 

operator Vij/wj is symmetric, sec Appendix. 

When the film contains several vortices positioned at 
some of the grid points r^, then the sum over these rj 
has to be performed in the first term of Eq. (28), and the 
total energy of this vortex system becomes: 

F = J2Fs{ri)+J2 Vij + E ('^') 

i j>i i 

* ^E^^^+E^-(^^)' (30) 

where the sums arc over the vortex positions z, j, Fs{ri) « 
^Vii is the self energy of the vortex, Vij is the vortex inter- 
action (29), and T4(ri) = J2i K^Ha{ri) is the potential 
caused by the applied field Ha{v). For constant Ha the 
external potential 14 (fj) = Ha Ylii has the shape of a 
negative trough which is zero along the film edges. The 
vortices are thus pulled into the film by this potential. 

G. Self energy of a vortex in the film 

In contrast to vortices in large bulk superconductors, 
the self energy of a vortex in a thin film depends on 
the film size and shape and on the vortex position even 
when it is far from the film edges. Calculating this vor- 
tex energy from the 2D current density in the film and 
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FIG. 4: The magnetic field Hz{x,y) in the plane of a thin 
square with hole and slit, Ha > 0, I — (case of flux fo- 
cussing, see Fig. 1 top right and Fig. 2 top) for ideal screening 
A = (top) and for A/a = 0.01 (bottom) (60 x 60 grid points, 
only half the square is shown). Note that even such small A 
strongly changes Hz{x,y) inside the superconductor, which 
penetrates much farther than A. The corresponding proflles 
Hz{x, 0) are show in Fig. 5. 



the 3D magnetic stray field outside the fihn would be 
a formidable task. Fortunately, a much simpler calcula- 
tion is possible using our above results and the known 
Lorentz force ^oVg{r) on the vortex: In a thought ex- 
periment, we move a first vortex from the film edge to 
its final position r^. At the edge the energy of this vor- 
tex, and also the interaction with other vortices needed 
later, are zero. Then its energy increases according to 
the integrated Lorentz force that originates from the in- 
teraction of this vortex with its own sheet current (more 
precisely: with the film edges, or with its images if an im- 
age method can be used, but this argument will not be 
required here) . When the vortex has reached position 
its energy is just its self energy Fs{ri). Now move a sec- 
ond vortex from the edge and merge it with the first one. 
The self energy of this new vortex, 4Fs{ri), is composed 
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FIG. 5: The profiles of the magnetic field Hz(x, 0) in the thin 
square of Fig. 1 (top right) taken along the x axis that passes 
through the hole and slit, in units of the applied field Ha 
and shown for several values of the 2D magnetic penetration 
depth A = \^/d = 0, 0.01, 0.03, 0.1, 0.3, and 1 in units of 
the half width a of the square. Same case as in Fig. 4 but 
with more grid points (100 x 100). Note that in the center of 
the square hole the magnetic field is enhanced by a factor of 
3 when A = 0, H^ (0, 0) ~ 3-ffa (flux focussing) ; in the narrow 
slit iifz(0.4a, 0) = 7AHa is even higher. Finite A reduces this 
enhancement and the spatial variation of H{x,0). 

of the two self energies and the energy needed to move 
the second vortex against the sheet current of the first 
one, equal to the interaction energy Vu, Eq. (29). From 
this we obtain the self energy Fs{ri) — 1^^/(4 — 2) = ^Vu 
used in Eq. (30). 

This result is exact within our numerical method, but 
in real films the self energy depends on the logarithm of 
the vortex core radius ~ ^, the coherence length. In our 
numerics ^ is effectively replaced by some cut-off length 
of the order of the grid spacing. If required, an improved 
consideration of the vortex core is possible if its radius 
exceeds the local grid spacing. The core shape may then 
be taken from the GL solution for infinite filmsiii If ^ 
is smaller than the grid spacing, the correct self energy 
is slightly larger than ^Vu, by a position-independent 
constant. 



H. The fluxoid in films 

The fluxoid <i>t inside a given closed path S inside the 
film is defined as the magnetic flux through this loop plus 
the 2D penetration depth A times the path integral of the 
sheet current, $f //io = Jg dx dy Hz{x, y)-\-K§ dS J(x, y). 
When g{x,y) and Hz{x,y) are given on a rectangular 
grid, the integration path is conveniently chosen along a 
closed rectangle which runs in the middle between the 
grid points (App. A). The components Jx and Jy are 
then obtained from the difference of the values of g(x, y) 
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FIG. 6: Thin superconducting square (half width a) with cen- 
tral square hole (half width ai) in applied field Ba = fioHa > 
and with no circulating current (/ = 0, case of flux fo- 
cussing). Plotted are the magnetic field B{0) — noHz{0,0) 
in the center of the hole (top) and the magnetic flux <& in- 
side the hole (bottom), for squares without slit (left, versus 
ai/a) and with a narrow radial slit (right, versus 1 — a-i/a), 
like in Fig. 1, for several values of the 2D penetration depths 
A/a = . . . oo. For A = and ai/a —* 0, both the minimum 
field and the average field in the hole diverge, 2B{Q)/Ba ~ 
$/4af_Ba « 0.80a/ai, see also the corresponding Figs. 15 and 
16 of Ref. 15 for rings. For small holes, a\ <^ a, the presence 
of a slit reduces this field enhancement. 



at neighboring grid points, while the flux is the sum of 
WkHzixk, Vk) over all points k inside this loop (App. A). 
The fluxoid obtained in this way from our solutions is 
indeed independent of the chosen path to within 4 to 
5 significant digits, confirming thus that the solutions 
g{x,y) and Hz{x,y) are accurate also for finite A and 
even when A <C a, see Figs. 7 and 8 below. Surprizingly, 
the obtained A-dependent g and Hz are quite accurate 
even when A is much smaller than the spacing of the 
numerical grid (typically > a/50), down to A/a = 10~^. 
Besides this, the solutions for A = are very accurate, 
with Hz — inside the film, see Fig. 5. 




0.4 , 0.6 

/ a 



FIG. 7: The fiuxoid <l>f trapped in the square hole of a thin 
square without radial slit in the flux-focussing case of Fig. 6. 
This figure is very similar to Fig. 14 of Ref. for rings. For 
small and large A nearly the same limiting curves are reached 
as for rings. The dots show the large-A limit for rings with 
radii ai and a: Aa\aBa/^t ~ 2aia ln(a/ai)/(a^ — a\)^ For 
A = this coincides with the $ in Fig. 6 (bottom left). 



III. EXAMPLE: FLUX FOCUSSING 

To illustrate how this method works and to check pre- 
vious approximations of a slitted ring by a full ringf2i4 

I discuss here the case of flux focussing in some de- 
tail. Consider a thin square film of half width a with 
a central square hole of half width ai (ai < \x\ < a, 

II < \y\ < a) without slit, or with a narrow radial slit 
of width A <^ a extending along the x axis from x ^ ai 
to X = a (Fig. 1). The current / circulating around 
the hole is forced to zero either "artificially" by putting 
the stream function g{x,y) = everywhere outside the 
superconducting film, i.e., also in the hole and slit, or, 
naturally, by cutting a slit that makes / = 0. A uniform 
magnetic field Ba = l^-oHa applied along z is screened 
inside the film (or partly screened if the 2D penetration 
depth is A = A^/d > for thickness d < A) by a sheet 
current J = — z x Vg = V x (zg) — {dg/dy, —dg/dx) 
that circulates clockwise near the edges and anticlock- 
wise near the hole, see Fig. 1 (top right). This screening 
current causes a magnetic field in the hole (and in the 
slit) that can be much higher than Ha, see Figs. 4 and 5. 

This field enhancement (or flux focussing) is plotted 
versus the relative hole size ai/a in Fig. 6. For squares 
with no slit our numerics yields for small holes with 
ai/a <C 1 and ideal screening (A = 0) for the central 
field B(0) = fioHziO, 0) « 0.405^ a/ai, and for the mag- 
netic flux in the hole $ w 0.80(4afi3a) a/ai, i-e., both 
values diverge as 1/ai when ai — > 0. This result is very 
similar to the flux focussing in circular rings depicted in 
Figs. 15 and 16 of Ref. 0. In particular, compared to a 
ring with radius a and hole radius ai, the limit A = 
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FIG. 8; The fluxoid <l>f trapped in the square hole of the 
square of Fig. 7, but now with a narrow radial slit. Since the 
integration path of the fluxoid must run inside the supercon- 
ductor, the path has to cross the slit by a narrow bridge, 
through which no current flows since 7 = in this flux- 
focussing case. The plot shows the fluxoid when this bridge 
is chosen at a; = ai (where the slit enters the hole, top) and 
at X = a (where the slit exits the square, bottom). For large 
holes these two <l?f are similar, but for small holes the larger 
integration path (bottom) yields a larger fluxoid since the in- 
tegration includes the not negligible magnetic flux inside the 
narrow slit. For A = the fluxoid of the smaller path (top) 
coincides with the flux $ of Fig. 6 (bottom right) . 



(where the flux $ equals the fluxoid <i>f) is almost iden- 
tical both for B{0)/Ba and for the trapped flux plot- 
ted as the effective area A^s = ^i/Ba referred to the 
hole area. Namely, for small holes at A = one has 
for squares AeQ/Aal ~ 0.80a/ai and for circular rings 
Aes/iral = (8/7r2)a/ai = 0.81a/ai (exact resul*MJ^). 

The corresponding results for the fluxoid <i>t in squares 
without slit are shown in Fig. 7. As expected, $f is inde- 
pendent of the integration path around the hole, which 
was chosen as any concentric square of half width be- 
tween oi and a. Again, these $f are very similar to those 
for circular rings shown in Fig. 14 of Ref. 0. For A = 




FIG. 9: The current stream lines in a thin superconductor 
square (|a;| < a, \y\ < a) with radial slit running at y = from 
a; = to a; = a, with a superconducting bridge ai x — a2 = 
0.35a. Shown is the flux-focusing case Ha > 0, 7 = 0, i.e., no 
current crosses the bridge. A/a — 0.01. The dots mark the 
numerical grid of 80 x 80 points used for this plot. 
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FIG. 10: The fluxoid <l>f and magnetic flux <1> inside a con- 
centric square of half width passing through the bridge (at 
a; = 02) of the slitted square of Fig. 9, plotted versus a2/a 
for several values A/a. For better presentation "I>f is shown 
4 times larger than "1>, and $ is plotted versus 1 — 02/0. For 
A = one has <l>f = $. For A/a » 1 one has $ = ia^Ba 
(dashed line, full penetration), and "l>f ~ [a2/ af''^ a? Ba. 



(where $f — $) they agree closely, as discussed above. 
But even for large A a and all hole sizes one has for 
ringsi Aeff — 2{o? — a\) / \n{a / ax) , which is also a good 
approximation for the square, see the dots in Fig. 7. 

Our 2D method allows us to check this approximation 
of slit-free flux focussing by considering squares or rings 
with slit. The results for the square with slit are depicted 
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in Fig. 6 [B{0) and $ plotted versus 1 — ai/a] and Fig. 8 
($f for two different integration paths). One can see 
that for large holes these realistic B{0), <I>, and $f are 
similar to those of slit-free squares. However, for small 
holes ai/a < 0.2, the field enhancement is considerably 
weakened by the presence of a slit, in particular for small 
A/a. For A = the enhancement factor does no longer 
diverge as 1 /ai but tends to saturate (or possibly diverges 
very weakly, as one over some logarithm, as is the case 
for finite A in the absence of a slit). This means that 
the slit changes the screening currents near small holes 
considerably and thus reduces the field enhancement. As 
can be seen from the curves in Fig. 6, the presence of a 
slit has qualitatively the same effect as an increased value 
of A. This finding applies even though the width A of 
our slit was very small, A/a ~ 1/500 and less. 

By the same token, the presence of a slit changes the 
fluxoid <I>f , Fig. 8. Moreover, the fluxoid now depends on 
the integration path. Since this path has to run inside the 
superconductor, one has to bridge the slit by a narrow 
superconducting bridge. The current through this bridge 
by definition is / = in this flux-focussing case. It turns 
out that the resulting <I>f depends on the position of this 
bridge. In Fig. 8 the two extreme cases are shown when 
this bridge is chosen at a; = ai (where the slit emerges 
from the hole) and at a; = a (where the slit exits the 
square). For large holes these two choices yield similar 
$f, but for small holes the large integration path yields 
a larger fluxoid than the small path. 

One reason for this difference is that the fluxoid for 
the large path includes the magnetic flux inside the 
slit. For A = one can show that is very large 

and almost does not decrease when the slit width A is 
decreased. From the simplified model of two long parallel 
strips (length I ^ 2a) with borders at y = ±a and y = 
±A/2 in perpendicular field Ba, one finds for narrow slits 
the trapped flujiiS 

= 7raZBa/ln(8a/A) for A < a . (31) 

For our squares with small hole we put the slit length 
I a and estimate the flux in the slit as = 
Tra^iJa/ ln(8a/A), which depends weakly on the slit 
width A. Comparing this with the flux in small holes 
from above, <i)*^°'° = (8/7r)ai aSa, we get the ratio 

(j)'^!'* TT^ a a 

$1^= 81n(8a/A) ^^^^ 

when A/a is of the order 1/500. This means, for small 
holes with relative width ai/a < 0.15, the magnetic 
flux even in a very narrow slit exceeds the flux in the 
hole. This flnding explains why for small holes our flux- 
focussing results for slit-free and slitted disks differ con- 
siderably while they agree for large holes. 

To check this further, we computed the magnetic flux 
$ and the fluxoid $f for a square with no hole, but with 
a narrow radial slit ranging from a; = to a; = a on the 
a: axis. We bridge this slit by a narrow superconduct- 
ing bridge centered at a; = a2, < a2 < a, see Fig. 9. 



The contour within which $ and $t are calculated is a 
concentric square passing through this bridge at a; = a2. 
We consider the case of flux- focussing {Ha > 0, / = 0), 
thus the current through the bridge is zero. In Fig. 10 
the resulting $ and $f are plotted versus a2/a for A/a = 
0, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, and oo. Note that 
the scales differ by a factor of 4. As expected, for A = 
(ideal screening) one has = $, and for A/a ^ 1 (full 
penetration) $ = Aa^Ba- Interestingly, For A/a <C 1, 
one has approximately $f « <i> « {ca2/a)4:a^Ba, where 
the constant c « 0.08 slightly depends on the slit width. 
This proportionality of the flux to the slit length within 
the square path, confirms that the magnetic fiux in the 
narrow slit is approximately proportional to its length, 
and that Eq. (31) (derived for a long double strip) is a 
good approximation for our radial slit in the square. 



IV. CONCLUSION 

In this paper I presented a method that allows to cal- 
culate the 2D distributions of the sheet current 3{x,y) 
and magnetic field component Hz{x,y) (and, of course, 
the full 3D magnetic field) for thin fiat superconductors 
of arbitrary shape. If the film thickness is d < A (the 
London depth), our method accounts for the 2D mag- 
netic penetration depth A = A^/d, which may have any 
value, < A < oo. The sheet current is expressed by the 
scalar potential (or stream function) g{x,y), for which we 
list 10 useful properties in Sec. II A. The statics and dy- 
namics of superconductors in the Meissner state, or with 
pinned and depinning vortices described by a complex or 
nonlinear resistivity, can be calculated. It is shown how 
this is generalized to multiply connected film shapes, e.g., 
squares or disks with a hole or closed slot, or with sev- 
eral such holes, slots, slits. For individual 2D vortices in 
the film we give their mutual interaction and self energy, 
which both depend on the specimen shape. The cou- 
pling of part of the magnetic fiux of a moving vortex into 
a hole or slit is expressed in terms of the stream function 
g{x,y), which is computed and depicted for some basic 
examples. 

If a slitted square or circular disk is used as a SQUID, 
the applied magnetic field and applied currents induce a 
signal across the weak link that bridges the slit, and the 
moving vortices cause a noisy signali The SQUID signal 
in principle can be calculated by our 2D method, see 
Ref. for a detailed theory of such SQUIDs and for the 
ID problem of a long double strip that models a linear 
SQUID. 

As a useful example for application of our 2D method 
we consider in Sec. Ill the phenomenon of flux focussing, 
which occurs when the total current circulating in a 
square or circular disk around a central hole is zero, 
/ = 0. We compare the "ideal case" when the disk has 
no slit (/ = can then be achieved by appropriate mag- 
netic history) with the realistic situation where a radial 
slit forces / = 0. We find good agreement for large holes. 
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but for squares with small central hole and radial slit, 
flux focussing is reduced by this slit. More applications 
will be published. 

Acknowledgments 

Helpful discussions with J. R. Clem, D. Koellc, G. P. 
Mikitik, Rinke Wijngaarden, and E. Zeldov are gratefully 
acknowledged. This work was supported in part by the 
German Israeli Research Grant Agreement (GIF) No G- 
705-50.14/01. 

APPENDIX A: NUMERICAL TRICKS 

For the computation of the stream function g{r) = 

g{x,y) with symmetry g(x,—y) = g{x,y), —a < x < a, 
< y < 6, a 2D grid of n points r/c = {xk,yk) with 
weights Wk, k = 1 . . .n, is chosen that covers the basic 
area 2a x b (upper half of the rectangular film 2a x 2b) . 
For simplicity we chose a rectangular grid Vij = (xi , yj ) 
with i = 1...2nx and j = 1...%, e.g., the equidis- 
tant grid Xi = —a + [i — ^)a/nx, yj = {j — ^)b/ny 
with constant weights w = 2ab/n, n — 2nx'ny. A bet- 
ter choice is a noncquidistant grid that is denser near 
the edges of the thin film. For example, for a rectangu- 
lar film with rectangular hole with borders at x = ±ai, 
y — ±&i and a narrow slit along y = Q (sec Fig. 1), a 
possible choice for the yj and weights Wyj (and corre- 
spondingly for the Xi and Wxi) is yj = b\f{vj), Wyj = 
hf'{vj)/nyi, Vj = {j - \)/nyi for j = l...nyi, and 
Vj =bi + (b- bi)f{vj), Wyj = {b- bi)f'{vj)/{ny - Uyi), 
'"j = (j - ^yi ~ 5)/ ("y - "yi) for i = %i -I- 1 . . . Uy, with 
any function f{v) defined in < u < 1 and having a 
derivative f'{v) = df/dv that at = and w = 1 is zero 
or small. We choose, e.g., 

f{v) = {Sv"^ -2v^)c+{l-c)v, 
f'[v) = &v{l-v)c+l-c, (Al) 

with 0<c<l;c = means equidistant yj with constant 

weights Wyj , and c = 1 (highest acciiracy) means that the 
distances yj+i—yj and weights Wyj vanish linearly at y = 
0, y = bi, and y = b. One can show that J2j '^yj^iVj) ~ 

(p{y)dy for any sufiiciently smooth function tp{y). 
For our 2D grid rij = {xi,yj) the weights are Wij = 
WxiWyj. This 2D counting of grid points is required for 



graphics and for computing derivatives, e.g., dg{x,y)/dx 
and d'^g{x,y)/dx'^. However, for computation of 2D in- 
tegrals or for matrix operations, a ID counting of the 
grid points is required: = = {xk,yk), Wij = Wk, 
k = i + 2nx{j — 1) = 1 . ■ .n, n = 2nx'ny, and any 2D func- 
tion like g{x, y) is represented as a vector gk — g{xk, yk)- 
The magnetic moment of the film, Sec. 2A, is then 
m = J g{r)d'^r « Y.k9kWk- 

The magnetic flux through a closed rectangular path 
running between the grid points (along x and y val- 
ues obtained from the above formulae for Xi, yi by us- 
ing half-integer values for i, j, e.g., j = 7/2) is then 
^ = Mo Yli'k ''^kHz{xk, yk) where the sum is over all points 
k inside this loop. 

Particular attention requires the computation of the 
Laplacian acting on g{x,y) [e.g., in Eq. (13)], and 
computed by multiplication by a matrix V^; [e.g., in 
Eq. (15)]. This operator should contain the information 
that g{x,y) = outside and on the outer edges of the 
rectangular film |a;| < a, \y\ < b, and that at y = one 
has dg{x,y)/dy = due to the symmetry g{x,—y) = 
g{x,y). A 2D method that in principle applies to any 
2D grid rfc = {xk, yk) computes V|, as the inverse of the 
Green fimction G'(r,r') satisfying V^G'(r,r') = S{r — r') 
and the conditions that G(r, r') = for r on the outer 
edge of the rectangle 2a x 2b and having even symme- 
try with respect to y. This G(r, r') may be expressed 
by an infinite sum, with alternating signs, of functions 
In ]r — r' — Rm„|/(47r), where the Rmn arc the vectors of 
a rectangular lattice with spacings 4a and 2b. The result- 
ing matrix indeed works, however, it is less accurate (and 
takes much more computation time) than the simple ID 
method of computing V^g = d'^g/dx^ + d'^g/dy'^ from 
g{xi,yj) and the values g{xi±i,yj±i) at the four neigh- 
boring points. With our noncquidistant grid we need 
for this the formula for f"{xi) = d'^f/dx^ at a; = Xi. 
Writing f{xi-i) = /_, f{xi) = fo, f{xi+i) = /+, 
hi = Xi — Xi-i > 0, /i2 = Xi+i — Xi > 0, one has 

/"(-o)-/-^-/o^ + /+/(^. (A2) 
hi + h2 hih2 hi+h2 

The boundary and symmetry conditions for g{x, y) allow 
to define the required values of g on the grid lines lying 
one grid spacing outside the basic area —a<x<a, 
< y < b: xo = -2a + xi, x„^+i = 2a - 
2/0 = y-i, Vny+i = 2b -yn^, g{xo,yj) = g{xn^+i,yj) = 
9{xi, yny+i) = 0, g{xi, yo) = g{xi, yi). 
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